Process chain approach to high-order perturbation calculus for quantum lattice 

models 
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A method based on Rayleigh-Schrodinger perturbation theory is developed that allows to obtain 
high-order series expansions for ground-state properties of quantum lattice models. The approach is 
capable of treating both lattice geometries of large spatial dimensionalities d and on-site degrees of 
freedom with large state space dimensionalities. It has recently been used to accurately compute the 
zero-temperature phase diagram of the Bose-Hubbard model on a hypercubic lattice, up to arbitrary 
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I. INTRODUCTION 

Although rather simple, quantum lattice gas or spin 
models are known to give rise to complex strongly cor- 
related many-body physics. Currently the interest in 
these paradigmatic models of condensed matter physics 
is stronlgy amplified by the fact that it just becomes ex- 
perimentally feasable to "engineer" many of them with 
ultracold atoms in optical lattice potentials^. Quan- 
tum lattice models are composed of small elementary 
sub-systems (sites) , a spin or the occupation of a single- 
particle state, arranged in a lattice. The geometry of the 
lattice is reflected in the fact that only neighboring sub- 
systems are coupled to each other directly. In general the 
Hamiltonian is of the form 



(1) 



where the on-site term hi acts on the local state space of 
lattice site % only, while the coupling term operates 
on both sites i and j with the sum Ylnj\ running over all 
pairs of neighboring sites (ij) (bonds). 

As the state-space dimensionality of the full system 
increases exponentially with the number of sites M, gen- 
erally the treatment of models like ([1]) is a hard prob- 
lem. However, sometimes there exist suitable approxima- 
tion schemes giving accurate results in certain regimes. 
One of them is given by high-order series expansions, ob- 
tained by automation of the usual Rayleigh-Schrodinger 
perturbation calculus. Such expansions have proven to 
be a useful tool for the investigation of zero-temperature 
properties^. In the simplest case, the coupling terms 
Vij are considered as perturbation. Starting from a 
product \g) = fj, \ n i) = \{ n i}) °f local eigenstates |rij) 
with hi\rii) = Ei\rii), physical quantities like expectation 
values or static susceptibilities (characterizing the state 
evolving from \g) adiabatically when the perturbation is 
switched on) can be expanded in high-order power se- 
ries with respect to a dimensionless coupling parameter. 
Also quantum phase transitions^, i.e. abrupt changes of 
the ground-state structure, occurring when a certain pa- 



rameter a passes a critical value a c , can be inferred from 
such expansions. 

A widely and successfully used algorithm for the ap- 
plication of perturbation calculus to such lattice mod- 
els is given by the connected cluster expansion^: As 
a first step one has to determine all sub-sets C (con- 
nected clusters) of mutually connected bonds (ij) pos- 
sessing not more than ^ max elements, with z^ max being the 
largest order of the expansion to be considered. More- 
over, all connected sub-clusters (and sub-sub-clulsters, 
. . . ) of each of these connected clusters have to be iden- 
tified. The second step consists of applying the perturba- 
tion calculus iteratively to each of the small sub-systems 
He = SieC hi + J2(ij)ec ^ij corresponding to the con- 
nected clusters C, with C" containing all sites connected 
to the bonds of C. Finally, these results can be combined 
to give the desired expansion. 

However, when trying to apply the connected clus- 
ter formalism to the Bose-Hubbard model^, describing 
bosonic particles with short-range interaction moving in 
a lattice of single-particle orbitals (sites), one encoun- 
teres two difficulties: (i) Finding all connected clusters as 
well as their sub-clusters is a demanding task for three- 
dimensional lattice geometries (and even more for spatial 
dimensionalities d > 3 that might be interesting in or- 
der to probe the convergence to the mean-field limit), 
(ii) The dimensionalities T>c of the connected cluster 
state spaces (that have to be considered after the parti- 
cle number conservation of the cluster Hamiltonians He 
has already been taken into account) grow extremely fast 
with respect to the averaged lattice site occupation (fill- 
ing) n. As a consequence, an iterative evaluation of the 
perturbation series up to 8th (10th) order is hindered 
already for the moderate filling n = 5 (n — 3) by the 
fact that it requires a representation of the cluster state 
spaces on the machine^. While difficulty (i) appears for 
every quantum lattice model of large spatial dimension- 
ality d > 3, difficulty (ii) is noticed for systems with the 
relevant local on-site state-space dimensionalities being 
too large. In the following, a method based on Kato's for- 
mulation of Rayleigh-Schrodinger perturbation calculus 9 



2 



will be described that circumvents both difficulties (i) 
and (ii). Here, perturbative corrections are obtained as 
sums over chains of processes acting in the "classical" 
space containing the unperturbed states but not their 
superpositions. These process chains, in turn, are gen- 
erated from paths through the <i-dimensional lattice. As 
an example, using this approach one is able to accurately 
compute the Bose-Hubbard phase diagram at any inte- 
ger filling n and for spatial dimensionalities d = 2, 3, and 
larger—. An implementation of the method presented 
here is straightforward and can be accomplished from 
scratch with reasonable effort. The present approach is 
found to be related to the ones described previously in 
Refs. [Illfl2lfl3l where contributions to the perturbation 
series are equally expressed in terms of sequences of pro- 
cesses. Differences to these approaches lie both in the 
way contributions are organized— and in the generation 
of diagrams from paths through the lattice. 

For clarity, the process chain method will be intro- 
duced in terms of the simple Bose-Hubbard model^. The 
generalization of the method to more involved or just dif- 
ferent lattice Hamiltonians (like those of Heisenberg-type 
spin models) is, however, straightforward; later on the 
properties of a general model amenable to the approach 
will be sketched. The Bose-Hubbard Hamiltonian reads 



H, 



BH 



E 



u 



n t (fii - 1) + (si - n)fii 



with S|, b { , and n, ; = b\b i denoting the bosonic creation, 
annihilation, and number operator for a single-particle 
orbital located at site i. The first line of Eq. @ includes 
an on-site interaction characterized by the energy cost U 
for each pair of particles occupying the same site. More- 
over, it assigns the local potential energy £j — /j, = —fa to 
particles sitting at site i, including an overall chemical po- 
tential fi introduced to control the total particle number 
N. The terms of the second line implement the kinetics, 
being exhausted by tunneling of particles between neigh- 
boring sites. Although rather simple, this model pro- 
vides a quantitative description of ultracold bosonic alkali 
atoms in optical lattice potentials^-. It shows quantum 
phase transitions between a gapless, compressible super- 
fluid phase with (quasi) long-range order present at large 
values of the ratio J/U and various gapped incompress- 
ible Mott-insulator phases (at sufficiently small J/U) 
with exponentially decaying correlations, each character- 
ized by an integer filling factor n ~ ^2i(hi) /M (depend- 
ing on the chemical potential fi/U)^&. 

This paper is organized as follows: Section |TT] is de- 
voted to the general perturbation expansion. It is briefly 
reviewed how expectation values and static susceptibili- 
ties can be obtained via the computation of energy cor- 
rections. Starting from Kato's formulation of the vth 
energy correction— then a way to reduce the number of 
terms that have to be taken into account to a minimum 



(2) 



is described. This is the first step that has to be ac- 
complished on the machine. Finally, the perturbative 
corrections are written as sums over process chains. This 
formulation will serve as a fruitful starting point for cus- 
tomizing perturbation calculus to quantum lattice mod- 
els in the way developed in section Unl The main concep- 
tual step of section IIIII consists in considering groups of 
operations (to be visualized by diagrams) such that each 
sequence of the operations contained in a group/diagram 
gives a process chain contributing to the desired pertur- 
bative correction. It is explained that the generation of 
the relevant diagrams can be put down to the rather sim- 
ple task of generating paths through the lattice (this al- 
lows one to address also large d), and it is described how 
the evaluation of diagrams can be performed. These are 
the two final steps to be implemented on the computer. 
In order to give a comprehensive presentation of the ap- 
proach, the computation of perturbative corrections is 
first discussed in the context of the Bose-Hubbard model 
@. Following along the lines given by this instructive 
example the application of the method to the large class 
of quantum lattice models described at the end of section 
IIIII should be straightforward. Section lTVl summarizes the 
basic steps of the approach, before section [V] closes with 
concluding remarks. 



II. GENERAL PERTURBATION EXPANSION 

A. Problem 

Consider a system described by the time-independent 
Hamiltonian 



H = H + XV 



(3) 



consisting of an unperturbed part Hq that is already di- 
agonalizcd, 



H = Y,E e \e)(e\, 



and a perturbation 



V = J2Ve',eW)(e\ 



(4) 



(5) 



multiplied by a dimensionless factor A finally to be set 
equal to 1. The standard non-degenerate Rayleigh- 
Schrodinger perturbation calculus^ gives a power law 
expansion 



En = E„ 



\ 2 E^ 



(6) 



for the energy of the state \G) evolving adiabatically from 
an eigenstate \g) of the unperturbed Hamiltonian (e.g. its 
ground state) when the perturbation is switched on. In 
the following the question of convergence will not be dis- 
cussed, but an algorithm for the computation of such 
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series expansions for many-body quantum lattice models 
is deviced. 

Usually, one is not interested in high-order perturba- 
tive corrections to the state of a many-body system, con- 
taining a vast amount of mostly unwanted information. 
One rather wishes to compute selected quantities like ex- 
pectation values (G\A\G). It is known, however, that 
the computation of ground state expectation values and 
static susceptibilities can be reduced to the evaluation of 
energy corrections. Introducing the extended Hamilto- 



H AB = H + XV + xA + yB, 



(7) 



with operators A and B, a perturbative treatment of 
V' = XV + xA + yB yields the expansion 



E, 



Gal 



E (v,m,k) x » x m y k 



(8) 



for the energy of the perturbed eigenstate \Gab) evolving 
from | <?). [In the case of non-hermitian operators A one 
can consider instead the hermitian ones A^ = ^(A+A") 
and = i(i-it) such that A = The 
low-order coefficients with respect to x and y then give 
series in powers of A for the expectation value 



(G|i|G)=^A^^ 1 »°), 



and the static susceptibilities 

XAB 

and 



XAA 



2,0) 



(9) 



(10) 



(11) 



describing the linear response of (A), when the full 
Hamiltonian H is perturbed by B or A, respec- 
tively. One can obtain these relations, e.g., by using 
the Hellmann-Feynman theorem -^(ip(z)\H(z)\ip(z)) — 

(i>(z)\[£H(z)]\il>(z)) [valid if \ip(z)} is a normalized 
eigenstate of H(z)\ or by interpreting them as first or- 
der energy corrections to the unperturbed problem given 
by the full Hamiltonian H . Hence, solely by using a for- 
malism for the evaluation of energy corrections, one can 
obtain series expansions for many important quantities 
characterizing the system. 



B. Minimal expression for the vt\\ order energy 
correction 

The vih energy correction appearing in Eq. (J6]) is given 
by Kato's closed expression^ (see also Ref. Il5h 

E g )= ^ ^{s a " +1 VS a » ■■■VS a3 VS a2 VS ai } (12) 
(f-i) 



with the sum runmn g over all combinations of the 

v+1 non-negative integers such that au = v—\, 

tr{-}^E e (e|-|e),and 



S a 



-\g){g\ 

J2 |e)(e| 



e^g (E g -E s )° 



for a = 
for a > 1 



(13) 



Since on the r.h.s. of Eq. (Ti"2"]) there are always at least 
two <Xi equal to zero, by cyclic permutation under the 
trace and using 



S a S a = 



-S° for a = and a' = 
for a = and a' ^ 

or a and a' = 
S a+a ' for a ^ and a' ^ 



(14) 



the energy correction can always be expressed as a sum 
over expectation values with respect to the unperturbed 
state \g). This gives 

E i V) = E G {ak} {g\VS^-W---S^VS^V\g) (15) 

(u-l) 

with 2fc=i a k = ^ — 1 according to the constraint of 
the sum (|12p and with (not uniquely determined) weight 
factors G| Qfc } taking into account how often each matrix 
element is generated with positive and negative prefactor 
during the elimination of the trace. Below the short hand 

(a t , . . . , a 2 , a x ) = ( 5 |V^f • • • S a >VS ai V\g), (16) 

for the matrix elements appearing in the sum (|15[) will be 
used, with () = (g|V|g) denoting the first order energy 
correction. 

An example for an expression like (|15|) is given by the 
formula^ (see also Ref. Il5j ) 



4-) 



(a„_i, . . . , a 2 , ai), 



(17) 



where all G{ afc } are either zero or one, as it is encoded in 
the set of constraints J2k=i a k — s with s = 1,2, ... (y — 
2), additional to the requirement Efc=i a k > v ~ 1- A 
further way to obtain an expression of type (|15p simi- 
lar to formula (|17p is to start with the matrix element 
(1,1,..., 1, 1), with all ak = 1, and to generate succes- 
sively further matrix elements to be considered by apply- 
ing the recursive scheme described in Ref. |17I . 

Many matrix elements (a u —i, . . . , aa, cci) that appear 
in the sums (|15p or (|17p give identical contributions: 
Writing explicitly S° — —\g)(g\, each matrix element 
(a„_i, . . . , a.2, Oi\) breaks up into elementary matrix ele- 
ments (EME) (/?£, . . . , /?2, Pi), containing strictly positive 
(i.e. non-zero) integers /3j only. Thus, e.g., one has 

(1,1,0,0,3,0,2,1) =-(l,l)()(3)(2,l) (18) 

= (1,1,0,0,3,0,1,2) =-(l,l)()(3)(l,2) 

= (1,1,0,0,2,1,0,3) =-(1, 1)()(2, 1)(3) 

= (1,1,0,0,1,2,0,3) =-(1, 1)()(1, 2)(3) 
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22 
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53 


12 
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119 


26 
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278 


47 


10 


627 


97 


11 


1,433 


180 


12 


3,216 


357 


13 


7,253 


668 


14 


16,169 


1,297 


15 


36,062 


2,427 


16 


79,876 


4,628 


17 


176,668 


8,637 


18 


388,910 


16,260 


19 


854,493 


30,188 


20 




56,252 


21 




103,848 


22 




191,873 


23 




352,204 


24 




646,061 



TABLE I: Minimum number K v of matrix elements that have 
to be considered for the vth. order energy correction. The 
number of contributing terms is significantly reduced to K' v 
if the first-order energy correction vanishes, (g\V\g) = () = 0. 
The data for ^ > 12 is taken from Ref. [Isl 

Here two basic operations were applied leaving the 
expression unchanged, the permutation of EMEs 
and the "reflection" of an EME, . . . , /3 2 , 0i) — > 
(/3i,/?2, •••!/%)• The latter is allowed since both V and 
the S a are hermitian. Hence, one has families of ma- 
trix elements (a v -i, ■ ■ ■ , ot2, a.\) with all members giv- 
ing the same contribution. In contrast, different EMEs 
(Pe, . . . , @2, 0i) generally give different contributions, if, 
by convention, one considers EMEs differing just by a 
"reflection" to be identical. Accordingly, each family of 
equally contributing matrix elements (a v —i, . . . , ct2, a{) 
is uniquely determined by a set of EMEs. In order to 
take into account only one representative of each family 
the energy correction can be rewritten as 

E t ] = E G W K-i.-> (19) 

(u-l) 

with a minimum number of non- vanishing weight factors 
G?™-i. Obtaining such a minimal expression for the or- 
ders to be considered is the first problem that has to 
be solved on a computer for an implementation of the 
method described here. Since so far one is dealing with 
the general perturbation expansion, this step has to be 
performed once only. 

A routine (Rl) for the generation of a minimal set 
of matrix elements contributing to the sum (|19p and 



their weights can be based on one of two alternative ap- 
proaches. The first one is to generate all matrix elements 
appearing in an expression like (|15p by starting cither 
from formula (fT2|) or (fTT|) , and then to identify mem- 
bers of the same matrix element family by decomposi- 
tion into EMEs. The second approach is to generate all 
EMEs and all combinations of them describing families of 
matrix elements (a„_i, . . . , «2, oil) appearing in the sum 
(fl~5f . The weight factor G™" , of a given family can then 

be obtained from a simple expression^. The minimum 
number K v of matrix elements to be considered in ^th 
order is listed in table |U It is drastically reduced, if it is 
known that the first order energy correction (<7|V"|<7) = () 
vanishes. Generally, given knowledge about the vanish- 
ing of corrections in certain orders v that are, say, even, 
odd or smaller than a value vq can be used to reduce 
the number of matrix elements that has to be taken into 
account. For example, if all corrections appearing in or- 
ders smaller than i/q are known to vanish, then there is 
at most one relevant matrix element appearing in order 
v . 

C. Energy corrections as sums over process chains 

One can interpret each matrix element 
(«!/_!, . . . , «2> cki) appearing in Eq. (| 19[) as a weighted 
sum J2{e t } over Paths \g) -> |ei) -> |e 2 ) ->•••—> 
|e„_i) — * \g) in a "classical" space containing the 
unperturbed states |e) but not their superpositions. All 
paths lead from \g) back to \g) via v — 1 intermediate 
states | eft): By plugging definitions ([5]) and (fl~3f into 
Eq. (fro|) . one obtains 

a 3 , ai ) = ( 2 °) 

{ e fcl 

*---w^v e2 , ei w^v ei , g 

with 

ee -S afi 6 e , g + (1 - <y a , )(l - S e , g )[E g - E e ]- a . (21) 

In that way a formulation involving huge quantum me- 
chanical state spaces is avoided. 

However, usually we have to deal with several per- 
turbing terms at once, V = V (1) + V {2) + ■■■ with 
y(m) _ J2 e , e V^ ™ |e')(e|, and we wish to keep track of 
them independently. Therefore, it is convenient to refor- 
mulate Eq. (|20p once more, namely as a sum over process 
chains P each of them being given by an ordered se- 

quence V^leiXsl- K&'xWH, • • • , V&l \a){e„-i\ 
of basic processes V^J \e') (e\ leading from \g) back to \g) 
in the classical space of unperturbed states introduced 
above. One has 

K_ l5 . . . , a 2 , Ql ) = v 9 { 2!M::; l} • ■ ■ ( 22 ) 
p 

^ " e 2 e 2 ,ei " ei e lt g ! 
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or, in combination with Eq. (|T9l 



E i v) = EE Qt-l 



. . j»(o2)y(m 2 )w(ai)y(mi) 
rv e 2 e 2 ,ei rf ei v ei,g 



(23) 



for the i/th order energy correction. 

The strategy for the computation of energy corrections 
proposed here is to generate all process chains P appear- 
ing in a given order n in a first step, and then to compute 
the contributions arising from each chain according to the 
different sets {ctk} possessing non-vanishing weight fac- 
tors in the general perturbation expansion in the 
form given by Eq. (fT^j) . In the next section it is shown 
how both steps can be accomplished efficiently for a lat- 
tice system with short-range coupling. 



(a) O 

(M) 






(e) 



o 



(f) 



III. LATTICE SYSTEM 
A. Bose-Hubbard problem 

In this section, the above formalism is applied to the 
Bose-Hubbard Hamiltonian Q on a hypercubic lattice 
geometry. All terms that are diagonal with respect to 
the lattice-site occupation numbers will be considered 
as unperturbed problem 



E 



V 

— ITU 

2 y 



1)_£ 



rii, 



(24) 



and the remaining tunneling terms as perturbation 



v = -±y 



b}b, 



(25) 



Energies have been expressed in units of the positive in- 
teraction parameter U, such that J/U is identified to be 
the dimensionless coupling parameter. A basis of unper- 
turbed states |e) is given by the lattice-site occupation- 
number states 



\{m}) = n 



(ft) 



(26) 



Let us assume that the state we want to investigate is the 
one evolving adiabatically from the unperturbed state 



\g) = \{n t =gi}) 



(27) 



when the perturbation is switched on. If \g) denotes the 
unperturbed ground state its occupation numbers gi min- 
imize [(<?i — l)/2 — [iij\J\gi and read 





h 

(h - 1) or h 



if Hi/U < 

if (h - 1) < m/U < h (28) 
if Hi/U = h-l>0, 



with non-negative integers h. As long as the marginal 
case of integer fXi/U is avoided, this state is protected by 
an energy gap. 



FIG. 1: Typical diagrams characterizing sets of tunneling 
operations between neighboring sites in the two-dimensional 
square lattice that appear in the energy correction of 2nd 
and 4th order. Tunneling operations are denoted by arrows, 
lattice sites by circles. Diagrams contributing to the pertur- 
bation expansion are those that can be interpreted as a single 
closed path, i.e. all except (d). 



In the following first a thorough discussion of the com- 
putation of (actual) energy corrections is given as an in- 
structive example. Then it will be shown that this ex- 
ample already contains everything one needs also for the 
computation of other quantities of interest, like single 
particle correlations (b\bj), number correlations (fiifij), 
or the static susceptibility x- -t for the annihilation and 

a 0' a 

creation operators of the condensate mode, a and a . 
The divergence of the latter indicates the quantum phase 
transition from a Mott-insulator to a superflui d 10 ! 19 ' 20 . 
The approach can easily be generalized to more compli- 
cated or just different lattice models. A general model 
that is amenable to the procedure described below is 
sketched at the end of this section. 



B. Corrections to the energy 

The perturbation V to be considered consists of tun- 
neling processes, i.e. the annihilation of one particle at a 
given site in combination with the creation of one parti- 
cle at a neighboring site. Denoting a tunneling operation 
by an arrow, one can visualize sets of tunneling oper- 
ations graphically by drawing diagrams. Obviously all 
process chains P starting and ending at the same (ar- 
bitrary) unperturbed state \g) must contain the same 
number of creation and annihilation processes at each 
site. Hence, those diagrams contributing to the energy 
corrections contain closed paths only. In Fig. [T] typical 
diagrams describing sets of tunneling operations in the 
two-dimensional square lattice are sketched. The num- 
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(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) 

FIG. 2: All topological diagrams contributing to the energy 
correction of a hypercubic lattice in the leading orders 2, 4 
and 6. 

ber of arrows corresponds to the power of J/U to which 
the diagram contributes. In principle one can obtain all 
process chains P appearing in the general formula (|23|) 
by generating all diagrams (sets of tunneling operations) 
contributing to a given order in a first step, and by order- 
ing the operations of each diagram in all possible ways in 
a second step. Hence, before energy corrections can be 
evaluated all contributing diagrams have to be generated. 

It has been noted that only diagrams containing closed 
paths of tunneling operations contribute, since only these 
have the same number of creation and annihilation oper- 
ations at each site. In correspondence with the connected 
cluster theorem^, one can also show that only connected 
diagrams give a non-vanishing contribution to the energy 
correction, i.e. those diagrams that cannot be divided 
into two or more sub-diagrams with no lattice site in com- 
mon or, in other words, that can be interpreted to consist 
of a single closed path only. For example, the energy cor- 
rections stemming from the different process chains that 
can be obtained from diagram (d) of Fig. [1] must add 
up to zero. The basic idea for a proof of this statement 
goes as follows^: A connected diagram would equally ap- 
pear in the perturbation expansion for a system differing 
from the one considered here by the fact that it consists 
of two completely independent sub-systems that are not 
coupled to each other by tunneling and with one part of 
the diagram lying in each of them. In that case, how- 
ever, the perturbed state evolving from an unperturbed 
product state will obviously be a product state with re- 
spect to both decoupled sub-systems, and its energy will 
be the sum of both sub-system energies. Thus, contribu- 
tions to the energy depending in a non-additive way on 
the properties of both sub-systems cannot occur. 

Unless one is not dealing with a rather small system, 
computing perturbative corrections to an extensive quan- 
tity like the energy will generally involve too many dia- 
grams to be accomplished in reasonable time. Neverthe- 
less, it is possible to compute these corrections for a ho- 
mogeneous system with /i, = fi (or a system that shows 
a different kind of translational symmetry). In that case 
topologically identical diagrams — like, e.g., (b) and (f) 
of Fig. Q] — will give identical contributions. In Fig. [2] 
all types of topologically different diagrams appearing in 
the leading orders^ 2, 4, and 6 are plotted. The multi- 
plicity Mt of a topological diagram T is a weight factor 
being defined as the number of ways it can be embedded 



into the given lattice geometry. Note that disconnected 
diagrams would give rise to unphysical multiplicities in- 
creasing with a power larger than one with the system 
size. In the following diagrams like those of Fig. [1] con- 
taining operations that are located in the lattice will be 
called geographical diagrams, while diagrams like those 
of Fig. [2] that are characterized by topology only will be 
called topological diagrams. 



C. Computing high-order energy corrections 

In order to obtain all topological diagrams contribut- 
ing to the energy correction and their multiplicities on 
a computer, one needs two basic routines, R2a and R2b, 
that will also serve to evaluate corrections to other expec- 
tation values than that of the energy. The first routine 
(R2a) computes all paths through the lattice via neigh- 
boring sites starting from a given site i to another site 
j containing v steps. By choosing i = j and associat- 
ing each step with a tunneling operation, in that way 
one obtains all sets of tunneling operations, i.e. all ge- 
ographical diagrams, contributing in order v. However, 
two corrections have to be taken into account. First, for 
closed loops a path visiting s different lattice sites could 
equally be associated to each of the s — 1 sites differ- 
ent from i. Hence, such a diagram must be weighted 
by a factor of s . Second, it might happen that differ- 
ent paths contain exactly the same tunneling operations 
(just in a different order). An example for that is given 
by the diagram shown in Fig.[T](b) that — starting from 
the site all arrows are connected to — might be obtained 
by first moving vertically and then horizontally or vice 
versa. If two paths include the same tunneling opera- 
tions, only one of them should be taken into account. For 
that purpose a further routine (R2b) is needed that iden- 
tifies paths containing the same tunneling operations. 

For a homogeneous Bravais lattice all sites are equal 
and it suffices to consider just a single site i; moreover 
only topologically distinct diagrams give distinct contri- 
butions. The relevant topological diagrams and their 
multiplicities can be obtained by collecting geographi- 
cal ones of the same topology. Note that this step serves 
only to reduce the number of diagrams that have to be 
evaluated. Hence, the algorithm used to probe the topo- 
logical equivalence of two geographical diagrams does not 
need to be perfect. A very simple way of identifying iden- 
tical topologies is to enumerate the sites appearing in a 
geographical diagram by a single index in the order they 
appear the first time in an associated path. Then rou- 
tine R2b can be used to compare the diagrams. For the 
closed loop diagrams contributing to the ground state en- 
ergy, with site i not being distinguished from the others 
ones appearing in the diagram, this approach has to be 
improved by probing enumerations starting at different 
sites. 

One advantage of the diagram generation via paths 
consist in that fact that it is easily implemented, even 
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for high spatial dimensionalities d. Assuming, for exam- 
ple, a hypercubic lattice, it is not difficult to design a 
routine R2a that in principle works for arbitrary d and 
also practically allows to consider values of d well above 
3 that might be interesting to study the convergence to- 
wards meanfield behavior. 

Once all diagrams of a given order have been obtained, 
their contribution to the energy (|23p can be evaluated 
by a last routine (R3). Each sequence of the opera- 
tions contained in a diagram corresponds to a differ- 
ent process chain P, with the permutation of two iden- 
tical operations — as they appear, e.g., in Fig. [5] (c) 
- not giving a new process chain. Thus, the dia- 
grams (a), (b), and (c) of Fig. [2] give rise to 2! = 2, 
4! = 24, and 4!/(2!) 2 = 6 different process chains, re- 
spectively. By applying the processes of a given chain 
one after the other to a small array of occupation num- 
bers {rii} initialized with = gi, giving a sequence 
{ n i°''}' { n l }> { n ^}i ■ ■ •> one can compute (i) the matrix 
elements V^] ek [being -{J/U)(nf\nf ] + 1)) 1/2 for 
tunneling from site i to site j], (ii) the unperturbed en- 
ergy differences E ek — E g of the intermediate states [given 

by £i {nt\n^-l)-g i {g } -l)-^/U){nf ) -g i )% and 
(iii) whether an intermediate state \ek) equals \g) {i.e. 
whether — gi for all i) or not. Afterwards one can 
choose only those matrix elements (a„_i, . . . , a,i, ai) ap- 
pearing in the minimal expression (|19p that have a/- = 
if |e/c) = Iff) and a.k ^ if |e&) ^ \g) for all intermediate 
states k — 1,2,..., n— 1. Only these give a non- vanishing 
contribution to the energy correction (|23[) that now can 
be evaluated with the help of Eq. (f2Tj) employing the 
energy differences E ek — E g . 

The scheme described in the preceding paragraph is 
not affected by the filling (the averaged particle num- 
ber per site), since no representation of a quantum- 
mechanical state space is needed. Moreover, while ap- 
plying the process chain to a set of occupation numbers 
those unperturbed basis states |{?ii}) that are relevant for 
a given order of perturbation calculus (and only those) 
are generated automatically. 



D. Corrections to expectation values and static 
susceptibilities 

Perturbative corrections to an expectation value (A) 
(or a static susceptibility \a g) m ^th order of the per- 
turbation V can be obtained by computing energy correc- 
tions for the combined perturbation V' = XV+xA+yB in 
first order in x (and y) and in ^th order in A, as expressed 
in Eq. ([9]) [and Eq. (flOf ] of section [TTJ Hence, in order to 
compute such quantities, one can proceed exactly as be- 
fore. Now just one process A e ^ e \e')(e\ associated to the 
operator A = £ e , A e i >e |e')(e| (and another one associ- 
ated to B) has to be included in each process chain, in 
addition to v tunneling processes stemming from V. The 




FIG. 3: Typical diagrams appearing in the perturbation ex- 
pansion of correlation functions (b]bj) in a 2D square lattice. 
The dashed arrow is associated to an operator b\bj. Again, 
all diagrams contain closed loops only, with the solid arrows 
describing a path from site i to site j, and disconnected dia- 
grams like (c) give zero-contribution. 

weight factors G^"} entering the general perturbation 
expansion (|23[) are those referring to energy corrections 
of order v + 1 for the computation of expectation values 
(or v + 2 for susceptibilities). 

For the computation of a correlation function (b\bj), 
the additional operation to be taken into account is the 
transfer of a particle from site j to site i described by the 
operator b\b-. Denoting such an operation by a dashed 
arrow, one can again use diagrams to describe sets of 
operations. Typical diagrams are shown in Fig. [3J As 
before, the requirement that any sequence of the oper- 
ations contained in a diagram must lead from a given 
unperturbed state \g) back to it (i.e. that the number of 
particle annihilations equals that of particle creations at 
every site) ensures that only closed loops (containing the 
dashed as well as solid arrows) contribute. Since, more- 
over, again disconnected diagrams like (c) don't need to 
be considered, the tunneling operations stemming from 
V can be interpreted as a single path leading from i to 
j. (One immediately sees that the perturbative treat- 
ment in the tunneling term discussed here is not sensi- 
tive to single-particle correlations between sites that are 
more than v steps between neighboring sites apart.) Also 
the generation of diagrams via the generation of paths 
through the lattice can be accomplished in a similar way 
as before by using routine R2a. Paths lead now from i to 
j and must not be corrected by the weight s , since no 
other starting point than the distinguished site i will be 
taken into account. Note that if i and j are neighboring 
sites, as it is the case in diagram (d) of Fig [31 the oper- 
ations related to the dashed arrow still gives a factor of 
x rather than A, such that it is well distinguished from a 
parallel tunneling operation stemming form V described 
by a solid arrow. Therefore, e.g. in diagram (d) the 
permutation of both upward tunneling operations (the 
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FIG. 4: Typical diagrams appearing in the perturbation ex- 
pansion of the number correlations (fuhj) in a 2D square 
lattice. The operation associated with the linked diamonds 
is described by the operator fctSj&j&j that leaves occupation 
numbers unaltered. 



"solid" and the "dashed" one) in a sequence does lead to 
a new process chain, such that here 4!/2! = 12 different 
process chains have to be taken into account. 

Another example is the computation of number cor- 
relations (hihj). The corresponding operator fiihj — 
bJbjbJS-, with matrix elements depending on the occu- 
pation of both sites i and j, leaves any unperturbed oc- 
cupation number state unaltered. Nonetheless one can 
associate this "operation" with a diagrammatic symbol 
that is chosen to be given by two diamonds at sites i 
and j, connected by a line. Fig. [4] shows some diagrams 
appearing in the perturbative expansion of expectation 
values (fiifij). Obviously, the tunneling operations must 
form closed paths, such that the number of creation and 
annihilation operations at each site are equal. Since, 
moreover, again only connected diagrams contribute, the 
closed tunneling paths must visit either i or j. For given 
i and j, one can generate all contributing diagrams by 
generating all combinations of two paths, such that one 
leads from i back to i and the other from j back to j 
(including paths of zero length). 

In the spirit of the examples treated so far, single-site 
expectation values (nf) with arbitrary power p are ob- 
tained from diagrams that are generated by paths leading 
from site i back to site i. 

It is worth mentioning that the expectation value of an 
operator like hifij, acting on a few sites only, can be com- 
puted even in the case of a large inhomogeneous system 
(provided, of course, perturbation theory is meaningful). 
Since the contributing geographic diagrams are just those 
exploring the neighborhood of i and j, their number is 
limited and does not depend on the system size. Correla- 
tions between i and j that are induced by the perturbing 
coupling term (like ■) in the present case) will, how- 
ever, only be taken into account in orders of perturbation 
theory that are comparable to the distance between i and 



FIG. 5: Diagrams contributing to the static susceptibility 
X- -t that indicates the quantum phase transition from a 
Mott-insulator to a superfluid. 




(x) (. ) { : ) (x) (x) 
® (*) » W 

(a) (b) (c) (d) (e) (f) (g) (h) 

FIG. 6: All topological diagrams contributing to the static 
susceptibility y„ -t in the leading orders 0, 1, 2 and 3. 



j (measured in steps between neighboring sites). This di- 
rectly reflects the limitation of the perturbative approach 
to systems with such correlations decaying on a distance 
being at most equal to the order of perturbation theory 
(at least as long as additional extrapolation techniques 
are not applicable or considered). 

Finally, it shall be outlined how the static susceptibil- 
ity x- -t = X f° r the annihilation and creation operators 

a ,a 

of the condensate mode, a and aj, with ao oc J2i f° r 
the homogeneous system can be computed. This quan- 
tity is proportional to the contribution oc |£| 2 to the 
energy obtained from the effective Hamiltonian H e g = 
H + + The process chains contributing 

to it contain (apart from tunneling processes) one cre- 
ation process and one annihilation process that will be 
represented diagramtically by a bullet • and a cross x , 
respectively. Examples for relevant diagrams are given 
in Fig. [5] They can be obtained from connected paths 
starting and ending anywhere in the lattice [including 
the zeroth order contribution shown in diagram (c)] . All 
topological diagrams appearing in the leading orders 0, 
1, 2, and 3 are shown in Fig. [51 In the limit of large spa- 
tial dimensionalities d, tunneling several times along the 
same bond becomes very unlikely such that diagrams like 
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(a), (b), (d), and (h) of Fig. [6] that can be interpreted as 
a path visiting each site only once give the major contri- 
bution to the susceptibility. 

For d > 2 the susceptibility x diverges when 
J/U reaches some critical value (J/U) c , indicating the 
quantum phase transition from a Mott-insulator to a 
superfluid 19,20 . However, the approximate value of \ 
in vt\\ order, x = 12k=o @ k (J/U) k with the coeffi- 
cients Pk depending on \i/U , will always be finite for 
finite J/U as long as /i/U is non-integer. Thus, in or- 
der to extract the critical parameter (J/U) c one has to 
resort to extrapolation to infinite order v. The critical 
parameter can be associated with the radius of conver- 
gence of the series for x with respect to J/U, namely 
(J/U) c = limfc^oo /3fe_i//3fc, assuming all /3k to have the 
same sign. Plotting ftk-i/flk versus 1/fc, one finds to 
good approximation all data points to lie on a straight 
line, suggesting the very simple phenomenological ex- 
trapolation scheme to extend the line to 1/fc — by 
a linear ffb21. This procedure gives the phase bound- 
ary (J/U) c versus \i/U with an estimated error of 1-2 
percent for arbitrary large filling n and spatial dimen- 
sionalities d = 2, 3, and greater^. The errors have 
been estimated by monitoring deviations of the approx- 
imate phase boundary while successively taking into ac- 
count more and more coefficients For n = 1 these 
results agree with those obtained by a strong coupling 
expansio n 2128 (d ~ 2) and by Quantum Monte Carlo 
simulation o 22 i 23 {d = 2 and 3). This simple example 
illustrates that extrapolation can be a valuable tool aug- 
menting high-order perturbation calculus. A brief intro- 
duction to more advanced extrapolation techniques as 
well as further references can be found in chapter f of 
Ref. H 

E. More general quantum lattice models 

So far, in this section the method has been developed 
in terms of the Bose- Hubbard model @. However, the 
approach is not restricted to this model, and following 
along the lines of the above example, it can be applied 
to a variety of quantum lattice models of the form given 
by Eq. ([I]). This includes fermionic Hubbard models or 
Heisenberg type spin models. In the remaining part of 
this section the properties of quantum lattice models that 
are amenable to the method described in this paper will 
be sketched. 

First of all the splitting of the full Hamiltonian ([I]) 
into an unperturbed part Hq and a perturbation V does 
not necessarily have to be such that Ho contains just 
all on-site terms hi, while V covers all coupling terms 
Vij. Both terms Hq and V can contribute to both on- 
site terms hi and coupling terms My. Moreover, the site 
index i can be generalized to run over several degrees 
of freedom at every lattice point, or, similarly, the sum 
J2{ij) can be extended to include not only pairs of near- 
est neighbors, but also further pairs of near sites. It is, 



however, required that the relevant set of eigenstates of 
H is characterized by a set of on-site quantum num- 
bers {m} taking values m = nf la , nf in + 1, rcf ln + 2, 
. . . , ax (with the possibility of arbitrary large on-site 
state-space dimensionalities T>i — n' nax — nf m + 1). In 
other words, Hq should be expressed in terms of number 
operators hi with hi\niri2 ■ • • Tii ■ ■ ■) = rii\nin2 ■ ■ ■ n.j • • •). 
For each site i there should also be a pair of ladder op- 
erators, if and l~ , being defined by if\nin 2 • • • n» • • •) = 
rff\n\n2 • ■ • ra» ± 1 • • ■). The perturbation V can be ex- 
pressed in terms of both number and ladder operators. 
Given the structure described in this paragraph, it will 
be possible to define diagrams and to evaluate them in a 
similar fashion as described for the Bose-Hubbard model 
above. 

Lattice systems covered by the scheme just outlined 
are bosonic or fermionic Hubbard models as well as spin 
Hamiltonians. In the former case the rii are just occu- 
pation numbers runing from n' nin = to n™ ax = oo for 
bosons and nf ax = 1 for fermions. The index i can also 
distinguish between different internal degrees of freedom 
(or species) of particles. For spins, n, would be associ- 
ated with the magnetic quantum number characterizing 
the spin at site i along a distinguished quantization axis, 
taking values between ±S with half-integer total spin S. 
Concerning just the implementation, the systems (or the 
unperturbed states) amenable to the approach described 
here do not need to be homogeneous and they can be de- 
fined on various lattice geometries; also frustration, dis- 
order, and certain types of long-range interaction can be 
present. But, of course, apart from being implementable 
the perturbation expansion must as well be a suitable 
approximation scheme for a given problem. 

For particles (i.e. in the Hubbard case) the unper- 
turbed Hamiltonian Hq can contain site-dependent po- 
tential terms cx hi and two-particle density-density inter- 
action terms cx hifij, also three- and more particle terms 
are possible. In spin models, corresponding terms can be 
considered, describing, e.g., local magnetic fields along 
the quantization axis or Ising type coupling. Non-local 
density-density interaction terms appearing in the un- 
perturbed Hamiltonian Hq can, in fact, be long-ranged. 
In that case the unperturbed energies computed during 
the evaluation of a given diagram will depend on the 
unerturbed quantum numbers m at sites not contained 
in that diagram (these will be unaltered by the opera- 
tions contained in the diagram). Considering a homoge- 
neous unperturbed state, this will cause only little ex- 
tra computational effort, for an inhomogeneous state the 
additional effort will just grow as the number of sites 
within the range of interaction. Clearly, the perturbative 
approach is limited to such strongly correlated phases 
that can be explored by starting from an unperturbed 
product state. However, the treatment of the bosonic 
Mott-transitio n 10 i 21 is an example showing that even the 
boundaries of such phases in parameter space can be ob- 
tained by applying suitable extrapolation schemes. One 
should also note that (near) degeneracies between the un- 
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perturbed state considered and other unperturbed states 
can spoil the results obtained by perturbation theory. 
As a simple example, this can happen for varying on-site 
potentials (or magnetic fields) that cause at some sites 
two different quantum numbers to lead to similar unper- 
turbed energies. As a remedy, in these cases one might 
consider to include degenerecy breaking terms in Hq and 
to subtract them again in the perturbation V, cf. chapter 
8 of Ref. 1 and references therein. 

The perturbation V can contain the ladder operators 
if. For the Hubbard models these correspond to the 
bosonic or fermionic creation and annihilation operators, 
for a spin model they are given by the raising and low- 
ering operators for the given quantization axis at site 
i. For bosonic and spin models the factors rff just de- 
pend on the local quantum number m. In the case of 
the fermionic Hubbard model, the factors T]f accompa- 
nying the creation or annihilation of particles take the 
values +1 and -1, depending on all occupation numbers 
{rii} (according to a given convention for the ordering 
of all sites i). Taking care of these signs will cause ad- 
ditional effort. While the unperturbed Hamiltonians Hq 
can contain long-range coupling terms, the coupling be- 
tween different sites i and j appearing in the perturbation 
V should be rather short-ranged, since the number of di- 
agrams to be evaluated grows rapidly with the number of 
coupling terms contained in V. If the coupling between 
different sites i and j is of the familiar form ifij, dia- 
grams can, again, be generated conveniently by finding 
paths through the lattice as described above for the Bose- 
Hubbard model. This form is, however, quite typical, as 
it describes both hopping of particles as well as spin-spin 
coupling in spin directions transverse to the quantization 
axis. 



IV. SUMMARY 

Let us briefly recapitulate the three basic steps of the 
approach described above in sections |TT] and IIIII The 
first task to be accomplished is to generate the leading 
order energy corrections (|19[) as they appear in standard 
Rayleigh-Schrodinger perturbation theory such that in 
every order v only a minimum number of different matrix 
elements (|16[) . each characterized by v — 1 non- negative 
integers ctk, has to be taken into account. This can be 
done, e.g., by starting from Kato's expression (|12p . and 
merging matrix elements that (via a decomposition into 
elementary matrix elements) are identified to give identi- 
cal contributions. The number of relevant terms can fur- 
ther be reduced if a priori knowledge is available about 
the vanishing of all matrix elements appearing in certain 
orders v with, e.g., v being even, odd, or smaller than 
some value v . The obtained results also serve for the 
computation of expectation values and static susceptibil- 
ities. 

The contributing matrix elements are interpreted as 
sums over process chains in a classical state space con- 



taining only the unperturbed states and not their super- 
positions, cf. Eq. (J22). For the lattice problems consid- 
ered, this formulation allows one to organize the pertur- 
bation expansion in terms of simple connected diagrams, 
each representing a collection of different operations. The 
second and the final third step to be accomplished are 
the generation and the evaluation of these diagrams. It 
has been shown that the generation of diagrams can be 
put down to the generation of paths through the lat- 
tice, a rather simple task that can be done on a com- 
puter even for large spatial dimensionalities d. Finally, 
the evaluation of diagrams is straightforward; one has to 
go through all possible sequences of the operations con- 
tained in a given diagram and map them according to 
Eq. (|22p to the terms of the general perturbation expan- 
sion (fTS)) obtained before. Since this procedure does not 
require a representation of the quantum- mechanical state 
space on the sub-lattice associated to a given diagram, 
it is not affected by large dimensionalities of the on-site 
state spaces. 



V. CONCLUSION 



A method to compute high-order series expansions for 
ground state properties of quantum lattice models has 
been described that is based on Rayleigh-Schrodinger 
perturbation calculus. The approach can be divided into 
three basic steps that have to be accomplished on a com- 
puter; each of them can be implemented with reasonable 
effort. Since the treatment of high spatial dimensionali- 
ties as well as of large lattice-site state-space dimensional- 
ities is not connected to serious difficulties, the presented 
approach complements the well-known connected cluster 
method^. Recently, the method described here has been 
used to compute the phase diagram of the Bose-Hubbard 
model on a d-dimensional hypercubic lattice, describing 
ultracold bosonic atoms in optical lattices^. It allowed 
not only to monitor in detail the convergence towards 
both the quantum-rotor limit of high filling n and the 
meanficld limit of large d, but also provided experimen- 
tally relevant data for two and three dimensional systems 
at moderate filling n = 2 — 10. However, as outlined 
in section IIIIEI a wide class of quantum lattice models, 
including Heisenberg-type spin and Hubbard-type tight- 
binding models are amenable to the approach described 
in this paper. These models can be frustrated and inho- 
mogeneous and they can contain disorder as well as long- 
range interaction of the density-density or Ising type. Es- 
pecially in view of the enormous interest in quantum lat- 
tice systems made of ultracold atomsi, the ease of treat- 
ing three-dimensional systems can make the method a 
valueable tool for current research. 
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In order to compute corrections to the ground state en- 
ergy in vih order, one has to consider connected clusters 
containing up to Mc = v sites, such that T>c is given by 
the number of possibilities to place Nc = nMc indistin- 
guishable particles on these sites. If about 500 bytes are 
needed per state space dimension, this gives 100 gigabyte 
of required RAM already for the filling of either n = 3 in 
10th order or n — 5 in 8th order. 

By taking Kato's expression (I12|l as a starting point here 
subcluster subtractions (or cumulant corrections) do not 
have to be performed on the level of each single diagram, 
as it is the case in Refs. [ll| and H3 . In contrast, the ap- 
proach briefly sketched in Ref. [HI is more similar to the one 
presented here; it does, however, not include the (huge) re- 
duction of terms contributing to the general perturbation 
expansion that will be described in subsection III Bl 
For any bipartite lattice geometry, as the hypercubic one, 
only an even number of tunneling operations leads back to 
the initial state such that only energy corrections of even 
power in J/U (not depending on the sign of the tunneling 
matrix element) appear. This fact can also be exploited to 
reduce the number of matrix elements considered in the 
general perturbation expansion (|23|l . 

Since the slope of the line is relatively low, the coeffi- 
cients Pk resemble approximately those of a geometric se- 
ries. Thus, also an exponential fit to the /3k already gives 
a good estimate to the phase boundary. 
Though using series expansions in powers of J/U, the ap- 
proach of Ref. [2J is different from the one described here. 
Instead of computing the susceptibility \ to obtain the 
phase boundary, the authors compare the energy of a state 
with integer filling to the energies of defect states having 
one additional particle or hole. Moreover, they use a con- 
nected cluster expansion to obtain energy corrections. 



